// snow_combination2.cpp : 为了支持第三季积雪业务重新修改合成策略
//

/*
FY-4A的亚洲（20-60N, 30-160E）区域产品和融合的中国区域产品还是用PolarSteregraph投影，IMS的欧亚产品（20-90N, 0-190E），谢谢。
另外，云能否用浅灰，others用深灰。

IMS、fy4和融合数据在日的阶段输出polar坐标数据，以后月、季、年合成数据和面积指数都基于这个polar坐标数据进行计算。


1-5
重新生成mask
扇形裁剪出图

2020-11-26 增加对第三级（亚洲高山区）的支持

2021-3-12 根据第三季需求重新修改合成策略和绘图方法（qgis），部署通过docker的方式。
月数据合成策略，从日数据合成，至少要有20天，否则不合成。
季数据从月数据合成，必须有三个月
年数据从月数据合成，采用自然年1-12月，必须有12个月
删除融合积雪的代码
每天写一个日志

2023-1-10
(1)process All Level3 扫描日期的bug 
(2)不用spdlog写日志，改用Airflow调用和记录日志

*/


#include <iostream>
#include "ogr_spatialref.h"//17
#include "gdal_priv.h"
#include <ctime>
#include <fstream>
#include "./sharedcodes-87/wftools.h"
#include "./sharedcodes-87/ajson5.h"
#include "./sharedcodes-87/wgt.h"
#include "./sharedcodes-87/wdb.h"
#include "combinecondition.h"
#include "spdlog/spdlog.h"
#include "spdlog/sinks/daily_file_sink.h"
#include "spdlog/sinks/stdout_color_sinks.h"
#include "productparam.h"
#include "wGdalRaster.h"

using namespace std;
using std::cout;

//运行模式 自动模式，补做模式，补做只绘图模式（覆盖老图）
#define RUNMODE_AUTO 0
#define RUNMODE_BUZUO 1
#define RUNMODE_BUZUO_PLOTONLY 2

//单日标准数据值的意义 3,6 for fy4
#define DAY_PIXVAL_OUTSIDE 0
#define DAY_PIXVAL_NODATA 3
#define DAY_PIXVAL_CLOUD 6
#define DAY_PIXVAL_WATER 7
#define DAY_PIXVAL_LAND 8
#define DAY_PIXVAL_ICE 9
#define DAY_PIXVAL_SNOW 10

//多期合成分类	IMS原始值
//范围外值	-9999
//缺测	-999
//常年陆地 -900 这个值仅在距平产品中有定义，用于区别距平为0的地区
//一次雪都没有且出现水的情况多于陆地和冰	-907
//一次雪都没有且出现陆地的情况多于水和冰	0（陆地颜色）
//一次雪都没有且出现冰的情况多于水和陆地	-909
//出现一次积雪及以上的像素	>0
#define MDAYS_PIXVAL_OUTSIDE -9999
#define MDAYS_PIXVAL_NODATA -999
#define MDAYS_PIXVAL_LAND_JP   -900
#define MDAYS_PIXVAL_WATER -907
#define MDAYS_PIXVAL_ICE -909

const int COMBMODE_DAY2MON = 0;//数据合成类型，日到月
const int COMBMODE_MON2SEAYEA = 1;//数据合成类型，月到季或者年



string g_gdaltranslate = "gdal_translate";
string g_gdalwarp = "gdalwarp";
string g_plotcmd = "" ;//绘图命令

double g_pixelarea = 16;//4km*4km 像素面积

int g_runMode =  RUNMODE_AUTO ;//auto 自动模式 , buzuo 补做模式 , buzuo_plotonly 补做仅仅绘图模式
string g_program = "" ;//程序路径，没啥用
string g_jsonfile = "" ;//输入json
int g_buzuoYmd0 = 0 ;//补做开始日期（包含）
int g_buzuoYmd1 = 0 ;//补做结束日期（包含）
string g_buzuoPid = "";//补做的产品id

//mysql operation
wDb g_db;//mysql数据库操作

vector<ProductParam> g_level3ppVec;//level3 月季年产品处理参数
vector<ProductParam> g_jupingppVec;//距平产品处理参数
ProductParam g_fy4dayPP ,  g_imsdayPP;//fy4日产品处理参数，ims日产品处理参数
GParams g_gpp ; //用于从json读取全局参数
QhcDataGroup g_qhc ;//气候场数据集



/// 积雪融合的参数，应该是不需要了
double g_A = 2 * 6371.228 / 25;
double g_B = 3.1415926 / 4;
double g_C = 3.1415926 / 180;

//将Fy4 03-06日产品和ims原始数据处理成标准日产品数据格式 包括polar坐标变换
void processFy4ImsDaily(ProductParam& pp, string intype);

//从文件名提取ymd
void extractFilesYmd(vector<string>& pathvec, vector<int>& ymdvec, int ymdloc)
{
	int n = pathvec.size();
	for (int i = 0; i < n; ++i)
	{
		string filepath = pathvec[i];
		string filename = wft_base_name(filepath);
		string ymdstr = filename.substr(ymdloc, 8);
		ymdvec.push_back(wft_str2int(ymdstr));
	}
}

//打印当前时间
void printCurrentTime() {
	//#include <ctime>
	time_t time0;
	time(&time0);
	std::cout << asctime(localtime(&time0)) << std::endl;
}

//积雪面积指数结构
struct AreaIndex
{
	inline AreaIndex() :maskid(0), sumval(0) {};
	int maskid;//区域ID
	double sumval;//面积指数值
};


//2021-4-13 将合成的具体文件记录在文本文件中
void write_combfiles( string path1 , vector<string> files )
{//outpath.cmbfiles.txt
  string txtpath = path1 + ".cmbfiles.txt" ;
  FILE* pf = fopen( txtpath.c_str() , "w" ) ;
  if( pf == 0 )
  {
      cout<<"Error: failed to create file "<<txtpath<<endl ;
  }else
  {
      for(size_t i = 0 ; i < files.size() ; ++ i){
          fprintf(pf,"%d  %s\n" , i+1 , files[i].c_str() ) ;
      }
      fclose(pf) ;
  }
}

//计算积雪面积指数
int computeAreaIndexManyDays(string maskfile, string snonwfile,
	int ndays, vector<AreaIndex>& result)
{
	double areaArr[255];
	short hasmask[255];
	memset(hasmask, 0, sizeof(short) * 255);
	memset(areaArr, 0, sizeof(double) * 255);
	GDALDataset* maskds = (GDALDataset*)GDALOpen(maskfile.c_str(), GA_ReadOnly);
	if (maskds == 0)
	{
		cout << "error : failed to open mask file " << maskfile << endl;
		return 10;
	}
	int x0 = maskds->GetRasterXSize();
	int y0 = maskds->GetRasterYSize();
	GDALDataset* snowds = (GDALDataset*)GDALOpen(snonwfile.c_str(), GA_ReadOnly);
	if (snowds == 0)
	{
		cout << "error : failed to open snow file " << snonwfile << endl;
		GDALClose(maskds);
		return 20;
	}
	int x1 = snowds->GetRasterXSize();
	int y1 = snowds->GetRasterYSize();

	if (x0 != x1 || y0 != y1)
	{
		cout << "error : snow and mask size are different ." << x0 << "," << y0 << ";" << x1 << "," << y1 << endl;
		GDALClose(maskds);
		GDALClose(snowds);
		return 30;
	}
	int bsize = x0 * y0;
	char* maskbuffer = new char[bsize];
	short* snowbuffer = new short[bsize];
	if (maskbuffer == 0 || snowbuffer == 0)
	{
		cout << "error : can not alloc memory." << endl;
		GDALClose(maskds);
		GDALClose(snowds);
		if (maskbuffer) delete[] maskbuffer;
		if (snowbuffer) delete[] snowbuffer;
		return 40;
	}
	maskds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, x0, y0, maskbuffer, x0, y0, GDT_Byte, 0, 0, 0);
	snowds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, x0, y0, snowbuffer, x0, y0, GDT_Int16, 0, 0, 0);
	GDALClose(maskds);
	GDALClose(snowds);

	if (ndays == 1)
	{
		for (int it = 0; it < bsize; ++it)
		{
			int maskid = (int)maskbuffer[it];
			hasmask[maskid] = 1;
			if (maskid > 0 && snowbuffer[it] == 10)
			{
				areaArr[maskid] += g_pixelarea;
			}
		}
	}
	else {
		for (int it = 0; it < bsize; ++it)
		{
			int maskid = (int)maskbuffer[it];
			hasmask[maskid] = 1;
			if (maskid > 0 && snowbuffer[it] > 0 )//bugfixed 2021-3-23
			{
				areaArr[maskid] += g_pixelarea * snowbuffer[it];
			}
		}
	}


	delete[] maskbuffer;
	delete[] snowbuffer;

	for (int i = 1; i < 255; ++i)
	{
		if (areaArr[i] > 0.001)
		{
			AreaIndex ai;
			ai.maskid = i;
			ai.sumval = areaArr[i] / ndays;
			result.push_back(ai);
		}
		else if (hasmask[i] > 0)
		{
			AreaIndex ai;
			ai.maskid = i;
			ai.sumval = 0;
			result.push_back(ai);
		}
	}
	return 100;
}




//保存积雪面积指数
string saveSnowAreaIndex(string pid, vector<AreaIndex>& vec, int ymd8, vector<string>& marr, vector<string>& pkarr )
{
    string sql ;
	for (int i = 0; i < vec.size(); ++i)
	{
        sql += g_db.insertUpdateXlSqlStr("tb_xulie_data",pid, wft_int2str(vec[i].maskid)
            ,wft_double2str(vec[i].sumval), wft_int2str(ymd8),
            marr , pkarr ) ;
	}
    return sql ;
}




//生成绘图的日期标题
string makePlotTitle( 
    string type, //day,mon,sea,yea
    int ymd 
    )
{
//Jan. 1 2020
//Jan. 2020
//Spring 2020
//2020
    int year = ymd/10000 ;
    int month = (ymd%10000)/100 ;
    int day = ymd%100 ;
//    static string monstrArr[] = { "January","February","March","April",
//        "May","June","July","August","September","October","November","December" };
    static string monstrArr[] = { "Jan.","Feb.","Mar.","Apr.",
        "May","Jun.","Jul.","Aug.","Sep.","Oct.","Nov.","Dec." };
	static string seastrArr[] = { "Spring","Summer","Autumn","Winter" };

    if( type=="day" ){
        return monstrArr[month-1] + " " + wft_int2str(day) + " " + wft_int2str(year) ;
    }else if( type=="mon"){
        return monstrArr[month-1] + " " + wft_int2str(year) ;
    }else if( type=="sea"){
        if(month==3) return seastrArr[0] + " " + wft_int2str(year) ;
        if(month==6) return seastrArr[1] + " " + wft_int2str(year) ;
        if(month==9) return seastrArr[2] + " " + wft_int2str(year) ;
        if(month==12) return seastrArr[3] + " " + wft_int2str(year) + "/" + wft_int2str(year+1) ;//update 2021-4-13
        return  string("unknown ") + wft_int2str(year) ;
    }else if( type=="yea"){
        return wft_int2str(year) ;
    }else{
        return wft_int2str(ymd) ;
    }
}

//绘图程序，使用qgis绘图
void plotqgis(
    string polarfile, //输入tif文件 
    string title,     //标题
	string temfile ,  //qgs模板文件
	string outfile    //输出文件
	 )
{
    string cmd1 = g_plotcmd + " " + temfile + " \"" + title + "\" " + polarfile + " " + outfile ;
    spdlog::info("try to call {}",cmd1) ;
    int ret = system( cmd1.c_str() ) ;
    spdlog::info("return val {}",ret) ;
}



//处理多天加和数据
//判断合成依据在前面做了，这里来了就只做合成，合成分两种，日到月，月到季或者年
//成功返回0
int computeSnowDaysAndOther(
    vector<string>& selfiles, //经过筛选后的文件数组，该函数内部不再进行筛选
	string outpath, //输出结果文件路径
    int combMode //合成方法，分两种，一个日到月 COMBMODE_DAY2MON ; 月到季或者年 COMBMODE_MON2SEAYEA
    )
{
	cout << "making sum " << outpath << " ... " << endl;
    spdlog::info("making {}" , outpath ) ;
	if (selfiles.size() == 0) {
        spdlog::warn("pick files empty.") ;
        return 2;
    }
	//get output parameters.
	int rasterXSize = 0;
	int rasterYSize = 0;
	{
		string dspath = selfiles[0];
		GDALDataset* tds = (GDALDataset*)GDALOpen(dspath.c_str(), GA_ReadOnly);
		if (tds == 0)
		{
			return 10;
		}
		rasterXSize = tds->GetRasterXSize();
		rasterYSize = tds->GetRasterYSize();
		GDALClose(tds);
	}
	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset* outputDataset = driver->Create(outpath.c_str(), rasterXSize, rasterYSize, 1, GDT_Int16, 0);
	if (outputDataset == 0)
	{
		cout << "Error : failed create " << outpath << endl;
        spdlog::warn("failed to create {}" , outpath) ;
		return 10;
	}
	short* buffer = new short[rasterXSize*rasterYSize];//输入文件的buffer，重复利用
	int bufferSize = rasterXSize * rasterYSize;
    vector<int> icecount(bufferSize,0) ;
    vector<int> landcount(bufferSize,0);
    vector<int> snowcount(bufferSize,0);  //累加积雪类型像素的次数（对于季和年不是累加天数，累加天数使用msnowcount）
    vector<int> msnowcount(bufferSize,0); //用于累加积雪天数，感觉好像也没啥必要分成两个变量，先留着吧（2021-3-17）
    vector<int> watercount(bufferSize,0);
    vector<int> outspacecount(bufferSize,0);
    vector<int> nodatacount(bufferSize,0);
    vector<int> availcount(bufferSize,0) ;//除了外部控件，其他类型总次数
    
	if (buffer == 0 )
	{
		spdlog::critical("processOneSum error : can not alloc more memories input buffer.");
		if (buffer) delete[] buffer;
		GDALClose(outputDataset);
		wft_remove_file(outpath);
		return 30;
	}

	int dsCount = selfiles.size();
	for (int ids = 0; ids < dsCount; ++ids)
	{
		string dspath = selfiles[ids];
		GDALDataset* dstemp = (GDALDataset*)GDALOpen(dspath.c_str(), GA_ReadOnly);
		if (dstemp == 0)
		{
            spdlog::critical("processOneSum error: failed to open '{}'.", dspath) ;
			std::cout << "*** Error : can not open " << dspath << ". out." << std::endl;
			return 31;
		}
		int nx = dstemp->GetRasterXSize();
		int ny = dstemp->GetRasterYSize();
		if (nx != rasterXSize || ny != rasterYSize) {
            spdlog::critical("processOneSum error: raster size different ({},{}),({},{}) '{}'.",
                rasterXSize, rasterYSize, nx, ny , dspath) ;
			std::cout << "*** Error : raster size is wrong " << rasterXSize << ","
				<< rasterYSize << ";" << nx << "," << ny << ". out." << std::endl;
			GDALClose(dstemp);
			return 32;
		}
		dstemp->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, rasterXSize, rasterYSize,
			buffer, rasterXSize, rasterYSize, GDT_Int16, 0, 0, 0);
		if (combMode == COMBMODE_MON2SEAYEA )
		{//month to season or year
			for (int it = 0; it < bufferSize; ++it)
			{
                if( buffer[it] == MDAYS_PIXVAL_OUTSIDE ) outspacecount[it]++ ;
                else if( buffer[it]==MDAYS_PIXVAL_NODATA ) nodatacount[it]++ ;
                else if( buffer[it]==MDAYS_PIXVAL_WATER ) watercount[it]++ ;
                else if( buffer[it]==MDAYS_PIXVAL_ICE ) icecount[it]++ ;
                else if( buffer[it]==0 ) landcount[it]++ ;
                else if( buffer[it]>0 ) {
                    snowcount[it]++ ;
                    msnowcount[it] += buffer[it] ;
                }
                if(buffer[it]!=MDAYS_PIXVAL_OUTSIDE) availcount[it]++ ;
			}
		}
		else
		{//day to month 
			for (int it = 0; it < bufferSize; ++it)
			{
                if( buffer[it] == DAY_PIXVAL_ICE ) icecount[it]++ ;
                else if( buffer[it]==DAY_PIXVAL_LAND) landcount[it]++ ;
                else if( buffer[it]==DAY_PIXVAL_OUTSIDE) outspacecount[it]++ ;
                else if( buffer[it]==DAY_PIXVAL_SNOW){
                    snowcount[it]++ ;
                    msnowcount[it]+=1 ;
                }
                else if( buffer[it]==DAY_PIXVAL_WATER) watercount[it]++ ;
                else if( buffer[it]==DAY_PIXVAL_CLOUD) nodatacount[it]++ ;
                else if( buffer[it]==DAY_PIXVAL_NODATA) nodatacount[it]++ ;
                else{//others unknow type use nodata
                    nodatacount[it]++ ;
                }
                if(buffer[it]!=MDAYS_PIXVAL_OUTSIDE) availcount[it]++ ;
			}
		}
		GDALClose(dstemp);
	}
   	vector<short> outBuffer(bufferSize) ;//输出文件的buffer
    for(int it = 0; it < bufferSize; ++it){
        if( outspacecount[it] > 0 ){//只要出现一次外部值就认为无效，设为外部值
            outBuffer[it] = MDAYS_PIXVAL_OUTSIDE ;
        }else if( nodatacount[it] > availcount[it]*0.2 ){
            outBuffer[it] = MDAYS_PIXVAL_NODATA ;//无效值缺测值个数大于总数20%，认为无效
        }else if( msnowcount[it] > 0 ){
            outBuffer[it] = msnowcount[it] ;
        }else if( watercount[it] > icecount[it] && watercount[it] > landcount[it] ){
            outBuffer[it] = MDAYS_PIXVAL_WATER ;
        }else if( icecount[it] > watercount[it] && icecount[it] > landcount[it] ){
            outBuffer[it] = MDAYS_PIXVAL_ICE ;
        }else {
            outBuffer[it] = 0 ; //land
        }
    }
	outputDataset->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, rasterXSize, rasterYSize,
		outBuffer.data(), rasterXSize, rasterYSize, GDT_Int16, 0, 0, 0);
	delete[] buffer; buffer = 0;
	{
		string dspath = selfiles[0];
		GDALDataset* tds = (GDALDataset*)GDALOpen(dspath.c_str(), GA_ReadOnly);
		double trans[6];
		tds->GetGeoTransform(trans);
		outputDataset->SetGeoTransform(trans);
		outputDataset->SetProjection(tds->GetProjectionRef());
		GDALClose(tds);
	}
	GDALClose(outputDataset);
	if (wft_test_file_exists(outpath) == false)
	{
		cout << "error : failed to make " << outpath << endl;
        spdlog::critical("processOneSum error: failed to write file '{}'",outpath);
		return 40;
	}else{
        return 0;//success
    }
}


//fy4数据赋予全盘坐标
void warpToGeos(string diskfile, string geosfile)
{
	string lonstr = "104.7";
	if (diskfile.find("995") != string::npos)
	{
		lonstr = "99.5";
	}
	//::gdal_translate -a_srs "+proj=geos +h=35785863 +a=6378137.0 +b=6356752.3 +lon_0=104.7 +no_defs" -a_ullr -5496000 5496000 5496000 -5496000 fy4.tif fy4.geos.tif
	{
		string cmd1 = g_gdaltranslate + " -a_srs \"+proj=geos +h=35785863 +a=6378137.0 +b=6356752.3 +lon_0=" + lonstr + " +no_defs\" "
			+ " -a_ullr -5496000 5496000 5496000 -5496000 " + diskfile + " " + geosfile;
		cout << cmd1 << endl;
		cout << system(cmd1.c_str()) << endl;
	}

	if (wft_test_file_exists(geosfile) == false)
	{
		cout << "failed to make " << geosfile << endl;
		return;
	}
}

//欧亚和亚洲使用95，中国使用105
void warpToPolar(string coorsysfile, string polarfile, string midlonstr = "95")
{
	//gdalwarp -t_srs "+proj=stere +lat_0=90 +lat_ts=0 +lon_0=95 +k=1 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs" -tr 4000 4000 -te -3800000 -4600000 3800000 1000000 ims.tif ims.pol.tif
	{
		string cmd1 = g_gdalwarp + " -overwrite -t_srs \"+proj=stere +lat_0=90 +lat_ts=0 +lon_0=" + midlonstr + " +k=1 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs\" "
			+ " -tr 4000 4000 -te -4050000 -4600000 4050000 1000000 " + coorsysfile + " " + polarfile;//modifyed
		cout << cmd1 << endl;
		cout << system(cmd1.c_str()) << endl;
	}
	if (wft_test_file_exists(polarfile) == false)
	{
		cout << "failed to make " << polarfile << endl;
		return;
	}
}





//将ims 0-4的值转换为统一格式0,7,8,9,10的值 
//增加判断全图是否有异常值，如果有类似205这样的异常值，跳过补做，检查trans中的分辨率，如果不是4000，跳过
//通过检查数据可用，返回true，否则返回false
bool convertIms5to10AndCheckExceptionValueAndResolution(string inpath, string outpath)
{
	GDALDataset* inds = (GDALDataset*)GDALOpen(inpath.c_str(), GA_ReadOnly);
	if (inds == 0)
	{
		cout << "failed to open " << inpath << endl;
        spdlog::critical("file:{} failed to open" , inpath) ;
		return false ;
	}
    double intrans[] = {0,0,0, 0,0,0};
    inds->GetGeoTransform( intrans ) ;
    if( fabsf(intrans[1]) != 4000 || fabsf(intrans[5]) != 4000 ){
        spdlog::critical("file:{} has invalid resolution {},{}, skip" , inpath, intrans[1] , intrans[5]) ;
        GDALClose(inds);
        return false ;
    }
    int xsize = inds->GetRasterXSize();
    int ysize = inds->GetRasterYSize();
    vector<unsigned char> inbuffer(xsize*ysize) ;
    inds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, xsize, ysize, inbuffer.data(), xsize, ysize, GDT_Byte, 0, 0, 0);

    
    int asize = xsize * ysize;
    vector<unsigned char> buffer(asize);
    bool hasExceptionValue = false ;
    int exceptionValue =0 ;
    for (int it = 0; it < asize; ++it)
    {
        //单日标准数据值的意义 3,6 for fy4
        switch (inbuffer[it])
        {
        case 0: buffer[it] = DAY_PIXVAL_OUTSIDE ;
        case 1: buffer[it] = DAY_PIXVAL_WATER; break;//sea
        case 2: buffer[it] = DAY_PIXVAL_LAND; break;//land
        case 3: buffer[it] = DAY_PIXVAL_ICE; break;//sea ice
        case 4: buffer[it] = DAY_PIXVAL_SNOW; break;//snow
        default:
            hasExceptionValue = true ;
            exceptionValue = inbuffer[it] ;
            buffer[it] = DAY_PIXVAL_NODATA;//new 2021-3-18 , 2004-11-12 appear 205 pixel value.
            break;
        }
    }
    
    if( hasExceptionValue==true ){
        spdlog::critical("file:{} has exception value {}, skip" , inpath, exceptionValue) ;
        GDALClose(inds);
        return false ;
    }
    
	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset* outds = driver->CreateCopy(outpath.c_str(), inds, 0, 0, 0, 0);
    GDALClose(inds);
    
	if (wft_test_file_exists(outpath)==true)
	{
		if (outds != 0)
		{
			outds->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, xsize, ysize, buffer.data(), xsize, ysize, GDT_Byte, 0, 0, 0);
			GDALClose(outds);
            return true ;
		}
		else {
			spdlog::critical("Error : failed to open {} ", outpath );
            return false ;
		}
	}
	else {
		spdlog::critical("Error : failed to make {} ", outpath );
        return false ;
	}
}


//将输入ims数据处理成标准polar7-10文件
void processImsPolar710( string inputfile, string outdayfile ,string outpolarfile ){
    if( g_runMode == RUNMODE_BUZUO_PLOTONLY ) return ;//补做绘图模式，不生产数据
    if (wft_test_file_exists(outdayfile) == false || g_runMode == RUNMODE_BUZUO )
    {
        //将IMS积雪产品值转换为5-10的统一值
        spdlog::info("begin processImsPolar710 {}", inputfile ) ;
        bool cok = convertIms5to10AndCheckExceptionValueAndResolution(inputfile , outdayfile);
        if( cok==true ){
            spdlog::info("processImsPolar710 {} good", inputfile ) ;
        }else{
            spdlog::critical("processImsPolar710 {} failed", inputfile ) ;
        }
    }
    //转统一polar投影
    if ( wft_test_file_exists(outdayfile) && 
        (wft_test_file_exists(outpolarfile) == false || g_runMode == RUNMODE_BUZUO) )
    {
        warpToPolar(outdayfile, outpolarfile);
    }
}
//将输入fy4数据处理成polar坐标
void processFy4Polar( string inputfile, string outpolarfile){
    if( g_runMode== RUNMODE_BUZUO_PLOTONLY) return ;//no production for plotonly
    
    string tempgeosfile = outpolarfile + ".geos.tmp.tif";
    warpToGeos( inputfile , tempgeosfile);
    if (wft_test_file_exists(tempgeosfile))
    {
        warpToPolar(tempgeosfile, outpolarfile);
        wft_remove_file(tempgeosfile);
    }else{
        spdlog::warn("failed to make:{}" ,tempgeosfile ) ;
    }
}



//确定产品是否需要补做
bool needBuzuo(const string buzuoid, const ProductParam& pp){
    if(  buzuoid == pp.xlpid )
    {
        return true ;
    }
    else{
        for(int ipl = 0 ; ipl < pp.plots.size() ; ++ ipl ){
            if( pp.plots[ipl].pid == buzuoid ){
                return true;
            }
        }
    }
    return false;
}

//process month level3 生产逐月产品
void processLevel3Month(const int scanyear0, const int scanyear1,
    const int scanymd0,const int scanymd1 ,const int biggestYmd,
    vector<int>& ymdVec , vector<string>& allfiles,
    ProductParam& pp
){
    for(int yearx = scanyear0 ; yearx <= scanyear1; ++ yearx ){
        for(int monx = 1; monx <= 12 ;++ monx ){
            int targetYmdi = yearx*10000+monx*100+1 ;
            if( targetYmdi < scanymd0 || targetYmdi> scanymd1 ) continue ;//skip
            //get ymd range 
            vector<string> pickfiles ;
            spdlog::info("try combine month {}-{}",yearx,monx) ;
            bool combok =isMonthCombineConditionOk(biggestYmd,ymdVec,allfiles,yearx,monx,pickfiles);
            if( combok==true ){
                for(int ii=0;ii<pickfiles.size();++ii) spdlog::info("{},{}",ii,pickfiles[ii]) ;
                
                string outname = pp.outprefix + wft_int2str(targetYmdi) + "." + pp.type + ".polar.tif";
                string outpath = wft_makeoutpath( pp.outdir , outname ,"yyyy" ) ;
                //2018-6-17 合成有几天合成几天，补做缺测判断，季和年必须满足3个月和12个月
                int oldcombnum = wft_read_combnumfile(outpath) ;
                //
                int newcombnum= pickfiles.size() ; 
                if(oldcombnum < newcombnum || g_runMode==RUNMODE_BUZUO || g_runMode==RUNMODE_BUZUO_PLOTONLY )
                {
                    if( g_runMode == RUNMODE_BUZUO_PLOTONLY){
                        spdlog::info("buzuo_plotonly mode , no data production for {}" , outpath) ;
                    }else{
                        spdlog::info("comb files count: old {} < new {} or in buzuo mode , need remaking {}.",
                            oldcombnum,newcombnum,outpath) ;
                        int statusSnow = computeSnowDaysAndOther(pickfiles, outpath, COMBMODE_DAY2MON);//当天像元值等于10时才加
                        if( statusSnow==0 ){
                            wft_write_combnumfile( outpath , newcombnum );
                            write_combfiles( outpath , pickfiles );//记录合成文件
                            //计算积雪面积指数 这里检查区域id为1的记录，如果没有就认为全部积雪指数都没有，要全部生产，另一个情况就是补做模式
                            if ( wft_test_file_exists(outpath) &&  (g_db.hasxl(pp.xlpid,"1",targetYmdi) == false 
                                    || g_runMode==RUNMODE_BUZUO) )
                            {
                                vector<string> maskidarr,pkarr ;
                                g_db.getxlmaskidandridarray(pp.xlpid,targetYmdi, maskidarr,pkarr ) ;
                                spdlog::info("get xl count {}, for xlpid {} ymd {}", maskidarr.size(), pp.xlpid, targetYmdi ) ;
                                vector<AreaIndex> aivec ;
                                int monthdaynum = cc_getMonthNDays(targetYmdi) ;
                                for(size_t ireg = 0 ; ireg < pp.maskfiles.size() ; ++ ireg  ){
                                    computeAreaIndexManyDays(pp.maskfiles[ireg],outpath,monthdaynum,aivec);
                                }
                                for(size_t imask=0;imask<aivec.size(); ++ imask){
                                    spdlog::info("ymd {},dtype {},mask {}:{}" , targetYmdi, pp.type, aivec[imask].maskid, aivec[imask].sumval) ;
                                }
                                string allsql = saveSnowAreaIndex(pp.xlpid, aivec, targetYmdi, maskidarr, pkarr );
                                int naffected = g_db.submitManySql(allsql) ;
                                if( naffected==-1 ){
                                    spdlog::critical("sql failed:{}" , allsql) ;
                                }
                                spdlog::info("naffected {}",naffected) ;
                            }
                        }else{
                            cout<<"error : computeSnowDaysAndOther retval "<<statusSnow<<endl ;
                            spdlog::critical("processLevel3Month error: computeSnowDaysAndOther return value {}",statusSnow);
                        }
                    }
                     
                    //区域绘图
                    for(size_t iplot = 0 ; iplot < pp.plots.size() ; ++ iplot ){ 
                        PlotParam& pl = pp.plots[iplot] ;
                        string pngfile1 = outpath.substr(0, outpath.length()-3 )+pl.tag+".png" ;
                        if (wft_test_file_exists(outpath) && 
                            (wft_test_file_exists(pngfile1)==false||g_runMode==RUNMODE_BUZUO||g_runMode==RUNMODE_BUZUO_PLOTONLY) )
                        {
                            string plotTitle = makePlotTitle( "mon" , targetYmdi) ;
                            plotqgis(outpath, plotTitle , pl.temfile ,pngfile1);
                        }
                        if(wft_test_file_exists(pngfile1) && 
                            (g_db.hasimg(pl.pid,targetYmdi)==false||g_runMode==RUNMODE_BUZUO||g_runMode==RUNMODE_BUZUO_PLOTONLY) )
                        {
                            g_db.insertUpdateImg( "tb_product_data" , 
                                pl.pid , outpath , wft_base_name(outpath) , 
                                pngfile1 , wft_int2str(targetYmdi) ) ;
                        }
                    }
                }
            }else{
                spdlog::info("combine condition not ok.") ;
            }
        }
    }
}


//process level3 season 生产季节产品 
void processLevel3Season(const int scanyear0,const int scanyear1,
    const int scanymd0 , const int scanymd1,const int biggestYmd,
    vector<int>& ymdVec , vector<string>& allfiles,
    ProductParam& pp
){
    for(int yearx = scanyear0 ; yearx <= scanyear1; ++ yearx ){
        for(int monx = 3; monx <= 12 ; monx+=3 ){
            int targetYmdi = yearx*10000+monx*100+1 ;
            if( targetYmdi < scanymd0 || targetYmdi> scanymd1 ) continue ;//skip
            //get ymd range 
            vector<string> pickfiles ;
            spdlog::info("try combine season {}-{}",yearx,monx) ;
            bool combok =isSeasonCombineConditionOk(biggestYmd,ymdVec,allfiles,yearx,monx,pickfiles);
            if( combok==true ){
                for(int ii=0;ii<pickfiles.size();++ii) spdlog::info("{},{}",ii,pickfiles[ii]) ;
                
                string outname = pp.outprefix + wft_int2str(targetYmdi) + "." + pp.type + ".polar.tif";
                string outpath = wft_makeoutpath( pp.outdir , outname ,"yyyy" ) ;
                if( wft_test_file_exists(outpath)==false || g_runMode==RUNMODE_BUZUO ||g_runMode==RUNMODE_BUZUO_PLOTONLY )
                {
                    if(  g_runMode==RUNMODE_BUZUO_PLOTONLY ){
                        spdlog::info("RUNMODE_BUZUO_PLOTONLY, no production for {}",outpath ) ;
                    }else{
                        spdlog::info("comb files count: {}, for {}",pickfiles.size(),outpath ) ;
                        int statusSnow = computeSnowDaysAndOther(pickfiles, outpath, COMBMODE_MON2SEAYEA);//
                        if( statusSnow==0 ){
                            //计算积雪面积指数 这里检查区域id为30的记录，如果没有就认为全部积雪指数都没有，要全部生产，另一个情况就是补做模式
                            if ( wft_test_file_exists(outpath) &&
                                (g_db.hasxl(pp.xlpid,"1",targetYmdi) == false || g_runMode==RUNMODE_BUZUO) )
                            {
                                vector<string> maskidarr,pkarr ;
                                g_db.getxlmaskidandridarray(pp.xlpid,targetYmdi, maskidarr,pkarr ) ;
                                spdlog::info("get xl count {}, for xlpid {} ymd {}", maskidarr.size(), pp.xlpid, targetYmdi ) ;
                                vector<AreaIndex> aivec ;
                                int seasondaycount = cc_getSeasonNDays( targetYmdi) ;
                                for(size_t ireg = 0 ; ireg < pp.maskfiles.size() ; ++ ireg  ){
                                    computeAreaIndexManyDays(pp.maskfiles[ireg],outpath,seasondaycount,aivec);
                                }
                                for(size_t imask=0;imask<aivec.size(); ++ imask){
                                    spdlog::info("ymd {},dtype {},mask {}:{}" , targetYmdi, pp.type, aivec[imask].maskid, aivec[imask].sumval) ;
                                }
                                string allsql = saveSnowAreaIndex(pp.xlpid, aivec, targetYmdi, maskidarr, pkarr );
                                int naffected = g_db.submitManySql(allsql) ;
                                if( naffected==-1 ){
                                    spdlog::critical("sql failed:{}" , allsql) ;
                                }
                                spdlog::info("naffected {}",naffected) ;
                            }
                        }else{
                            cout<<"error : computeSnowDaysAndOther retval "<<statusSnow<<endl ;
                            spdlog::critical("processLevel3Season error: computeSnowDaysAndOther return value {}",statusSnow);
                        }
                    }
                    
                    //区域绘图
                    for(size_t iplot = 0 ; iplot < pp.plots.size() ; ++ iplot ){ 
                        PlotParam& pl = pp.plots[iplot] ;
                        string pngfile1 = outpath.substr(0, outpath.length()-3 )+pl.tag+".png" ;
                        if (wft_test_file_exists(outpath) && 
                            (wft_test_file_exists(pngfile1)== false||g_runMode==RUNMODE_BUZUO||g_runMode==RUNMODE_BUZUO_PLOTONLY) )
                        {
                            string plotTitle = makePlotTitle( "sea" , targetYmdi) ;
                            plotqgis(outpath, plotTitle , pl.temfile ,pngfile1);
                        }
                        if ( wft_test_file_exists(pngfile1) && 
                            (g_db.hasimg(pl.pid,targetYmdi)==false||g_runMode==RUNMODE_BUZUO||g_runMode==RUNMODE_BUZUO_PLOTONLY) )
                        {
                            g_db.insertUpdateImg( "tb_product_data" , 
                                pl.pid , outpath , wft_base_name(outpath) , 
                                pngfile1 , wft_int2str(targetYmdi) ) ;
                        }
                    }
                    
                    
                }
            }else{
                spdlog::info("combine condition not ok.") ;
            }
        }
    }
}


//process level3 year 生产年产品
void processLevel3Year(const int scanyear0,const int scanyear1,
    const int scanymd0 , const int scanymd1,const int biggestYmd,
    vector<int>& ymdVec , vector<string>& allfiles,
    ProductParam& pp
){
    for(int yearx = scanyear0 ; yearx <= scanyear1; ++ yearx ){                
        int targetYmdi = yearx*10000+101 ;
        if( targetYmdi < scanymd0 || targetYmdi> scanymd1 ) continue ;//skip
        //get ymd range 
        vector<string> pickfiles ;
        spdlog::info("try combine year {}",yearx ) ;
        bool combok =isYearCombineConditionOk(biggestYmd,ymdVec,allfiles,yearx,pickfiles);
        if( combok==true ){
            for(int ii=0;ii<pickfiles.size();++ii) spdlog::info("{},{}",ii,pickfiles[ii]) ;
            
            string outname = pp.outprefix + wft_int2str(targetYmdi) + "." + pp.type + ".polar.tif";
            string outpath = wft_makeoutpath( pp.outdir , outname ,"yyyy" ) ;
            if( wft_test_file_exists(outpath)==false || g_runMode==RUNMODE_BUZUO||g_runMode==RUNMODE_BUZUO_PLOTONLY )
            {
                if(g_runMode==RUNMODE_BUZUO_PLOTONLY){
                    spdlog::info("RUNMODE_BUZUO_PLOTONLY, no production for {}",outpath ) ;
                }else{
                    spdlog::info("comb files count: {}, for {}",pickfiles.size(),outpath ) ;
                    int statusSnow = computeSnowDaysAndOther(pickfiles, outpath, COMBMODE_MON2SEAYEA);//
                    if( statusSnow==0 ){
                        //计算积雪面积指数 这里检查区域id为30的记录，如果没有就认为全部积雪指数都没有，要全部生产，另一个情况就是补做模式
                        if ( wft_test_file_exists(outpath) &&(g_db.hasxl(pp.xlpid,"1",targetYmdi) == false || g_runMode==RUNMODE_BUZUO) )
                        {
                            vector<string> maskidarr,pkarr ;
                            g_db.getxlmaskidandridarray(pp.xlpid,targetYmdi, maskidarr,pkarr ) ;
                            spdlog::info("get xl count {}, for xlpid {} ymd {}", maskidarr.size(), pp.xlpid, targetYmdi ) ;
                            vector<AreaIndex> aivec ;
                            int yeardaycount = cc_getYearNDays(targetYmdi) ;
                            for(size_t ireg = 0 ; ireg < pp.maskfiles.size() ; ++ ireg  ){
                                computeAreaIndexManyDays(pp.maskfiles[ireg],outpath,yeardaycount,aivec);
                            }
                            for(size_t imask=0;imask<aivec.size(); ++ imask){
                                spdlog::info("ymd {},dtype {},mask {}:{}" , targetYmdi, pp.type, aivec[imask].maskid, aivec[imask].sumval) ;
                            }
                            string allsql = saveSnowAreaIndex(pp.xlpid, aivec, targetYmdi, maskidarr, pkarr );
                            int naffected = g_db.submitManySql(allsql) ;
                            if( naffected==-1 ){
                                spdlog::critical("sql failed:{}" , allsql) ;
                            }
                            spdlog::info("naffected {}",naffected) ;
                        }
                    }else{
                        cout<<"error : computeSnowDaysAndOther retval "<<statusSnow<<endl ;
                        spdlog::critical("processLevel3Year error: computeSnowDaysAndOther return value {}",statusSnow);
                    }
                }
                //区域绘图
                for(size_t iplot = 0 ; iplot < pp.plots.size() ; ++ iplot ){ 
                    PlotParam& pl = pp.plots[iplot] ;
                    string pngfile1 = outpath.substr(0, outpath.length()-3 )+pl.tag+".png" ;
                    if (wft_test_file_exists(outpath) && 
                        (wft_test_file_exists(pngfile1)==false||g_runMode==RUNMODE_BUZUO||g_runMode==RUNMODE_BUZUO_PLOTONLY) )
                    {
                        string plotTitle = makePlotTitle( "yea" , targetYmdi) ;
                        plotqgis(outpath, plotTitle , pl.temfile ,pngfile1);
                    }
                    if ( wft_test_file_exists(pngfile1) && 
                        (g_db.hasimg(pl.pid,targetYmdi)==false||g_runMode==RUNMODE_BUZUO||g_runMode==RUNMODE_BUZUO_PLOTONLY) )
                    {
                        g_db.insertUpdateImg( "tb_product_data" , 
                            pl.pid , outpath , wft_base_name(outpath) , 
                            pngfile1 , wft_int2str(targetYmdi) ) ;
                    }
                }
            }
        }else{
            spdlog::info("combine condition not ok.") ;
        }
    }
}


//读取日数据的polar目录下文件，保存到level3的polar目录下，用以解决995和104.5不能一起统计的问题
///处理三级数据，就是月、季、年数据，输入数据为天数据的polar投影数据。统一使用polarstere投影，
/// 这样就不需要考虑995或者104.5或者其他的投影数据了。 
/// pole3ok
void processLevel3( vector<ProductParam>& ppvec )
{
	int todayYmd = wft_current_dateymd_int();
	for (int ipp = 0; ipp < ppvec.size(); ++ipp)
	{
		ProductParam pp = ppvec[ipp];
		cout<<"begin "<<pp.pname <<"-"<< pp.type <<endl ;
        spdlog::info("begin pname {}, type {}, enable {}" , pp.pname , pp.type ,pp.enabled) ;
		if( g_runMode == RUNMODE_AUTO )
		{
			if (pp.enabled != 1) {
				cout << pp.pname << " not enabled." << endl;
                spdlog::info("{} not enabled" , pp.pname) ;
				continue ;//skip
			}
		}else
		{
			if( needBuzuo(g_buzuoPid , pp) == false )
			{
                spdlog::info("{} not need buzuo", pp.pname) ;
				continue ;//skip
			}
		}
		
		vector<string> allfiles;
		vector<int> ymdVec;
        int todayYmd = wft_current_dateymd_int() ;
        int scanymd0 = todayYmd; 
        int scanymd1 = todayYmd;
        //对接的数据目录是3-6点合成的fy4日产品
        if( g_runMode == RUNMODE_AUTO){
            int todayYear = todayYmd / 10000 ;
            //2023-1-10 bug // scanymd0 = (todayYear-1)*10000+0101 ;
            scanymd0 = (todayYear-1)*10000+ 101 ;//bugfixed 
            scanymd1 = todayYmd ;
        }else{
            //补做模式
            scanymd0 = g_buzuoYmd0 ;
            scanymd1 = g_buzuoYmd1 ;
        }
        cout<<"scan for "<<scanymd0<<" - "<<scanymd1<<endl ;
        spdlog::info("runmode {} , scan ymd {}-{}" ,g_runMode, scanymd0, scanymd1);
        int scanyear0 = scanymd0/10000 ;
        int scanyear1 = scanymd1/10000 ;
        for(int iy = scanyear0 ; iy <= scanyear1; ++ iy ){
            string dir2 = pp.indir + wft_int2str(iy) ;
            wft_get_allSelectedFiles( dir2 , pp.inprefix ,  pp.intail , -1 , "" , allfiles ) ;
        }
		int numfiles = allfiles.size();
        cout<<"all scan files:"<<numfiles<<endl ;
        spdlog::info("all scan files {}", numfiles) ;
		if (numfiles == 0)
		{
			cout << "Warning : no files is need to process." << endl;
			continue;
		}
		wft_descSortStrVec( allfiles) ;
        
        int biggestYmd = 0 ;
		for (int i = 0; i < numfiles; ++i)
		{
			string filepath = allfiles[i];
			string filename = wft_base_name(filepath);
			string ymdstr = filename.substr(pp.ymdloc, 8);
			int ymdi = (int)atof(ymdstr.c_str());
			ymdVec.push_back(ymdi);
			if( ymdi > biggestYmd ) biggestYmd = ymdi ;
		}
        spdlog::info("biggest file ymd {}" , biggestYmd) ;
        
        if( pp.type=="mon" ){
            processLevel3Month( scanyear0 , scanyear1 , scanymd0 , scanymd1 , biggestYmd , 
                ymdVec , allfiles , pp ) ;
        }else if( pp.type=="sea"){
            processLevel3Season( scanyear0 , scanyear1 , scanymd0 , scanymd1 , biggestYmd , 
                ymdVec , allfiles , pp) ;
        }else if( pp.type=="yea"){
            processLevel3Year( scanyear0 , scanyear1 , scanymd0 , scanymd1 , biggestYmd , 
                ymdVec , allfiles , pp) ;
        }else{
            spdlog::info("unknown pp.type '{}'." , pp.type) ;
        }
	}
}


//处理FY4每日积雪
//intype: fy4/ims
void processFy4ImsDaily(ProductParam& pp, string intype)
{
    spdlog::info("in processFy4ImsDaily") ;
	vector<string> allfiles;
	vector<int> ymdVec;
	string inprefix = pp.inprefix;
	string intail = pp.intail;
	int todayYmd = wft_current_dateymd_int() ;
    
    int scanymd0 = todayYmd; 
    int scanymd1 = todayYmd;
    
	
    //对接的数据目录是3-6点合成的fy4日产品
    if( g_runMode == RUNMODE_AUTO){
        int todayYear = todayYmd / 10000 ;
        int todayMonth = (todayYmd % 10000)/100 ;
        
        scanymd0 = todayMonth-1 ;
        if( scanymd0 == 0 ) scanymd0 = (todayYear-1)*10000+1201 ;
        else scanymd0 = todayYear*10000 + (todayMonth-1)*100 + 1 ;
        scanymd1 = todayYmd ;
    }else{
        //补做模式
        scanymd0 = g_buzuoYmd0 ;
        scanymd1 = g_buzuoYmd1 ;

    }
    cout<<"scan for "<<scanymd0<<" - "<<scanymd1<<endl ;
    spdlog::info("runmode {}, scan ymd {}-{}" ,g_runMode, scanymd0, scanymd1);
    int scanyear0 = scanymd0/10000 ;
    int scanyear1 = scanymd1/10000 ;
    for(int iy = scanyear0 ; iy <= scanyear1; ++ iy ){
        string dir2 = pp.indir + wft_int2str(iy) ;
        wft_get_allSelectedFiles( dir2 , inprefix ,  intail , -1 , "" , allfiles ) ;
    }

	int numfiles = allfiles.size();
	if (numfiles == 0)
	{
		cout << "Warning : no files is need to process." << endl;
        spdlog::info("no scan files.") ;
		return;
	}

	cout<<"sorting "<<allfiles.size()<<" files..."<<endl ;
	wft_descSortStrVec(allfiles) ;


	for (size_t i = 0; i < numfiles; ++i)
	{
		string filepath = allfiles[i];
		string filename = wft_base_name(filepath);
		string ymdstr = filename.substr(pp.ymdloc, pp.ymdlen);
        if( pp.ymdlen==7 ){
            string yearstr = ymdstr.substr(0, 4);
            string doystr = ymdstr.substr(4, 3);
            int tyear = wft_str2int(yearstr);
            int tdoy = wft_str2int(doystr);
            int tmon, tday;
            wft_convertdayofyear2monthday(tyear, tdoy, tmon, tday);
            int tymd8 = tyear * 10000 + tmon * 100 + tday;
            ymdstr = wft_int2str(tymd8);
        }
		int ymdi = (int)atof(ymdstr.c_str());
		ymdVec.push_back(ymdi);
	}
    int per0 = -1 ;
	for (size_t i = 0; i < numfiles; ++i)
	{
		int ymdi = ymdVec[i];
        
        //补做或者自动模式筛选数据
        if( ymdi < scanymd0 || ymdi > scanymd1 ) continue;//skip

		//圆盘转极投影
		string outname2 = pp.outprefix + wft_int2str(ymdi) + ".day.polar.tif";
		string polarfile = wft_makeoutpath( pp.outdir ,  outname2 , "yyyy") ; 
        
        //补做绘图模式不生产数据
        if( g_runMode != RUNMODE_BUZUO_PLOTONLY ){
            if( intype=="fy4" ){
                processFy4Polar( allfiles[i] , polarfile) ;
            }else if( intype=="ims" ){
                string outdayname2 = pp.outprefix + wft_int2str(ymdi) + ".day.tif";
                string outdayfile = wft_makeoutpath( pp.outdir ,  outdayname2  , "yyyy") ; 
                processImsPolar710( allfiles[i] , outdayfile , polarfile ) ;
            }
            //计算积雪面积指数 这里检查区域id为30的记录，如果没有就认为全部积雪指数都没有，要全部生产，另一个情况就是补做模式
            if ( wft_test_file_exists(polarfile) &&  (g_db.hasxl(pp.xlpid,"1",ymdi) == false || g_runMode==RUNMODE_BUZUO) )
            {
                vector<string> maskidarr,pkarr ;
                g_db.getxlmaskidandridarray(pp.xlpid,ymdi, maskidarr,pkarr ) ;
                spdlog::info("get xl count {}, for xlpid {} ymd {}", maskidarr.size(), pp.xlpid, ymdi ) ;

                vector<AreaIndex> aivec ;
                //2021-3-12 这里使用全国的一个mask文件取代使用分省求和
                for(size_t ireg = 0 ; ireg < pp.maskfiles.size() ; ++ ireg  ){
                    computeAreaIndexManyDays(
                        pp.maskfiles[ireg] ,//mask文件
                        polarfile ,//输入文件
                        1 ,//几天
                        aivec);
                }
                for(size_t imask=0;imask<aivec.size(); ++ imask){
                    spdlog::info("ymd {},dtype {},mask {}:{}" , ymdi, pp.type, aivec[imask].maskid, aivec[imask].sumval) ;
                }
                string allsql = saveSnowAreaIndex(pp.xlpid, aivec, ymdi, maskidarr, pkarr );
                int naffected = g_db.submitManySql(allsql) ;
                if( naffected==-1 ){
                    spdlog::critical("sql failed:{}" , allsql) ;
                }
                spdlog::info("naffected {}",naffected) ;
            }
        }else{
            spdlog::info("RUNMODE_BUZUO_PLOTONLY no production for {}" , polarfile ) ;
        }

		//区域绘图
        for(size_t iplot = 0 ; iplot < pp.plots.size() ; ++ iplot ){ 
            //绘图
            PlotParam& pl = pp.plots[iplot] ;
            string pngfile1 = polarfile.substr(0, polarfile.length()-3 )+pl.tag+".png" ;
            if (wft_test_file_exists(polarfile) && 
                (wft_test_file_exists(pngfile1)==false||g_runMode==RUNMODE_BUZUO||g_runMode==RUNMODE_BUZUO_PLOTONLY) )
            {
                string plotTitle = makePlotTitle( "day" , ymdi) ;
                //plot image by qgis
                plotqgis(polarfile, plotTitle , pl.temfile ,pngfile1);
            }

            //绘图入库
            if ( wft_test_file_exists(pngfile1) && 
                (g_db.hasimg(pl.pid,ymdi)==false||g_runMode==RUNMODE_BUZUO||g_runMode==RUNMODE_BUZUO_PLOTONLY) )
            {
                g_db.insertUpdateImg( "tb_product_data" , 
                    pl.pid , polarfile , wft_base_name(polarfile) , 
                    pngfile1 , wft_int2str(ymdi) ) ;
            }
        }
        int per1 = i*100/numfiles ;
        if( per0!=per1 ){
            per0 = per1 ;
            cout<<per0<<" " ;
        }
	}
    cout<<"100."<<endl ;
}

//生产一个剧评文件，并绘图入库
void process1Juping(
    ProductParam& pp,
    const string infile,//输入数据文件
    const int ymdi,     //年月日 yyyyMMdd
    const string dtype, //时间类型 mon/sea/yea
    QhcDataGroup& qhc ,  //气候场数据集对象 qhc 
    string outjpfile    //输出剧评文件路径
){
    if( wft_test_file_exists(outjpfile) ==false || g_runMode==RUNMODE_BUZUO )
    {//没有文件或者补做模式，那么生产剧评产品tif
        wGdalRasterI16 indata ;
        bool inok = indata.open(infile) ;
        if( inok==false ){
            spdlog::warn("process1Juping, open infile failed {}.", infile ) ;
            return ;
        }
        
        int xsize = indata.getXSize() ;
        int ysize = indata.getYSize() ;
        int asize = xsize*ysize ;
        vector<short> qhcdata ;
        bool qhcok = qhc.getQhcData( dtype , ymdi ,qhcdata ) ;
        if( qhcok==false ){
            spdlog::warn("process1Juping, open qhcfile failed for type {} ymd {} .",dtype, ymdi ) ;
            return ;
        }
        
        if( asize != qhcdata.size() ){
            spdlog::warn("process1Juping, indata size {}  not equal qhcdata size {} .", asize , qhcdata.size() ) ;
            return ;
        }
        
        wGdalRasterI16 outdata ;
        bool outok = outdata.create(xsize,ysize,1) ;
        if( outok==false ){
            spdlog::warn("process1Juping, create outdata failed .") ;
            return ;
        }
        //geotrans and proj
        outdata.copyProj( indata.getProj() ) ;
        outdata.copyTrans( indata.getTrans() ) ;
        
        //逐像素处理
        for(int it = 0 ; it <asize;++ it ){
            int v = indata.getValuei(it,0) ;
            int q = qhcdata[it] ;
            if( q==0 && v==0 ){
                //常年陆地
                outdata.setValuei(it , 0 , MDAYS_PIXVAL_LAND_JP ) ;
            }else if( (q>0 && v>=0) || (q>=0 && v>0) ){
                outdata.setValuei(it , 0 , v-q ) ;
            }else if( q>=0 ){
                outdata.setValuei(it , 0 , v) ;
            }else if( v>=0 ){
                outdata.setValuei(it , 0 , q) ;
            }else{
                outdata.setValuei(it , 0 , q) ;
            }
        }
        
        //saving
        bool saveok = outdata.save( outjpfile ) ;
        if( saveok == false ){
            spdlog::warn("process1Juping, failed to save {}", outjpfile) ;
        }else{
            spdlog::info("process1Juping, save ok: {}", outjpfile) ;
        }
    }
    //区域绘图
    for(size_t iplot = 0 ; iplot < pp.plots.size() ; ++ iplot ){ 
        PlotParam& pl = pp.plots[iplot] ;
        string pngfile1 = outjpfile.substr(0, outjpfile.length()-3 )+pl.tag+".png" ;
        if (wft_test_file_exists(outjpfile) && 
            (wft_test_file_exists(pngfile1) == false || g_runMode==RUNMODE_BUZUO||g_runMode==RUNMODE_BUZUO_PLOTONLY) )
        {
            string plotTitle = makePlotTitle( pp.type , ymdi) ;
            plotqgis(outjpfile, plotTitle , pl.temfile ,pngfile1);
        }
        if ( wft_test_file_exists(pngfile1) && 
            (g_db.hasimg(pl.pid,ymdi)==false||g_runMode==RUNMODE_BUZUO||g_runMode==RUNMODE_BUZUO_PLOTONLY) )
        {
            g_db.insertUpdateImg( "tb_product_data" , 
                pl.pid , outjpfile , wft_base_name(outjpfile) , 
                pngfile1 , wft_int2str(ymdi) ) ;
        }
    }
}


//处理一个剧评产品
void processJupingProduct( ProductParam& pp, QhcDataGroup& qhc )
{
    spdlog::info("start processing {}" , pp.pname) ;   
    vector<string> allfiles;
    int todayYmd = wft_current_dateymd_int() ;
    int scanymd0 = todayYmd; 
    int scanymd1 = todayYmd;
    if( g_runMode == RUNMODE_AUTO){
        int todayYear = todayYmd / 10000 ;
        //2023-1-10 bug // scanymd0 = (todayYear-1)*10000+0101 ;
        scanymd0 = (todayYear-1)*10000+ 101 ;//fixed
        scanymd1 = todayYmd ;
    }else{
        //补做模式
        scanymd0 = g_buzuoYmd0 ;
        scanymd1 = g_buzuoYmd1 ;
    }
    cout<<"scan for " << scanymd0<<" - "<<scanymd1<<endl ;
    spdlog::info("runmode {}, scan ymd {}-{}" ,g_runMode, scanymd0, scanymd1);
    int scanyear0 = scanymd0/10000 ;
    int scanyear1 = scanymd1/10000 ;
    for(int iy = scanyear0 ; iy <= scanyear1; ++ iy ){
        string dir2 = pp.indir + wft_int2str(iy) ;
        wft_get_allSelectedFiles( dir2 , pp.inprefix ,  pp.intail , -1 , "" , allfiles ) ;
    }
    int numfiles = allfiles.size();
    cout<<"all scan files:"<<numfiles<<endl ;
    spdlog::info("all scan files {}", numfiles) ;
    if (numfiles == 0)
    {
        cout << "Warning : no files is need to process." << endl;
        return ;
    }
    wft_descSortStrVec( allfiles) ;
    for (int i = 0; i < numfiles; ++i)
    {
        string filepath = allfiles[i];
        string filename = wft_base_name(filepath);
        string ymdstr = filename.substr(pp.ymdloc, 8);
        int ymdi = (int)atof(ymdstr.c_str());
        if( ymdi < scanymd0 || ymdi > scanymd1 ){
            continue ;//skip
        }
        string outname = pp.outprefix + wft_int2str(ymdi) + "." + pp.type + ".polar.tif";
        string outjpfile = wft_makeoutpath( pp.outdir , outname ,"yyyy" ) ;
        process1Juping( pp, filepath , ymdi , pp.type , g_qhc , outjpfile) ;
        cout<<i<<"/"<<numfiles<<" " ;
    }
    cout<<"all done. "<<endl ;
}



int main(int argc, char** argv)
{ 
	printCurrentTime();
	std::cout << "Description: snow combination for  daily(day),month(mon),season(sea),year(yea). " << std::endl;
	std::cout << "Version 1a . by wangfengdev@163.com 2017-12-16." << std::endl;
	std::cout << "Version 2a . add snow fusion v3. by wangfengdev@163.com 2017-12-18." << std::endl;
	std::cout << "Version 2.0.1a . bugfixed 20171219-0757 ." << std::endl;
	std::cout << "Version 2.1.1a . do or not do  ." << std::endl;
	std::cout << "Version 2.2a . add mjzs for season and year. 2017-12-21  ." << std::endl;
	std::cout << "Version 3a . use ronghe v4 2017-12-31  ." << std::endl;
	std::cout << "Version 4.0 . big modify for production procedure, 2018-1-4  ." << std::endl;
	std::cout << "Version 4.1 . Ronghe Snow plot in albers proj map, 2018-1-16  ." << std::endl;
	cout<<"Version 4.2 buzuo for fy4 ,rnh ,ims daily product. 2018-3-31"<<endl ;
	cout<<"Version 4.2.1 add overwrite, modify level3. 2018-4-2"<<endl ;
	cout<<"Version 4.2.2 history mode. 2018-4-24"<<endl ;
	cout<<"Version 4.2.3 add fiv,ten. 2018-4-24"<<endl ;
	cout<<"Version 4.2.4 add fiv,ten. 2018-4-25"<<endl ;
	cout<<"Version 4.2.5 fusion. 2018-4-25"<<endl ;
	cout<<"Version 4.2.6 cbtics. 2018-4-26"<<endl ;
	cout<<"Version 4.2.7 fy4 daily. 2018-5-11"<<endl ;
	cout<<"Version 4.3 comb any day count. 2018-6-17"<<endl ;
	cout<<"v5.0.2 add third pole."<<endl ;
	cout<<"v6.0.2 big update for combination method, add anomaly product, map render using qgis, add daily log.2021-3-12"<<endl ;
    cout<<"v6.1.1 add log path, plot title bug fixed."<<endl ;
    cout<<"v6.2.9 generate and submit a group sql 2021-3-13."<<endl ;
    cout<<"v6.3  add db connection fail msg 2021-3-15"<<endl ;
    cout<<"v6.4.1 update combine method , juping method , SAI method. add buzuo_plotonly mode.  2021-3-17"<<endl ;
    cout<<"v6.4.2 ice water land compare bug fixed. 2021-3-17"<<endl ;
    cout<<"v6.5.1 process ims bad value, skip ims full image outspace, resolution 1km skip. 2021-3-18"<<endl ;
    cout<<"v6.5.2 jp product add proj and trans 2021-3-18."<<endl ;
    cout<<"v6.5.4 fix findxl(...) bad param of 30, fix AreaIndex bug for many days.2021-3-23"<<endl; //2021-3-23
    cout<<"v6.5.5 params=2 unknown mode bug fixed.2021-3-23."<<endl ;
    cout<<"v6.6  记录日到月具体合成的文件路径,修改冬季的标题2020/2021. 2021-4-13"<<endl ;
    cout<<"v6.6.1 scan bugfixed; no need spdlog instead use airflow. 2023-1-10"<<endl ;
	if (argc == 1)
	{
		std::cout << "*** sample call: ***" << std::endl;
		std::cout << "auto mode:          snow_combination2 startup.json" << endl;
		std::cout << "buzuo mode:         snow_combination2 startup.json buzuoPid buzuoYmd" << endl;
        std::cout << "buzuo mode:         snow_combination2 startup.json buzuoPid buzuoYmd0 buzuoYmd1" << endl;
        std::cout << "buzuo_plotonly mode:snow_combination2 startup.json buzuoPid buzuoYmd0 buzuoYmd1 plotonly" << endl;
        cout << "no enough params." << endl;
		return 11;
	}

	std::cout << "load params..." << std::endl;
	g_program = argv[0] ;
	g_jsonfile = argv[1];
	g_runMode = RUNMODE_AUTO;
	if( argc==4 )
	{//buzuo
		g_runMode = RUNMODE_BUZUO ;
		g_buzuoPid = argv[2] ;
		g_buzuoYmd0 = atof(argv[3]) ;
        g_buzuoYmd1 = g_buzuoYmd0 ;
        spdlog::info("start for buzuo {}." ,g_buzuoYmd0 ) ;
	}else if( argc==5 ){
        g_runMode = RUNMODE_BUZUO;
		g_buzuoPid = argv[2] ;
		g_buzuoYmd0 = atof(argv[3]) ;
        g_buzuoYmd1 = atof(argv[4]) ;
        spdlog::info("start for buzuo from {} to {}.",g_buzuoYmd0,g_buzuoYmd1) ;
    }else if( argc==6 ){
        g_runMode = RUNMODE_BUZUO_PLOTONLY;
        g_buzuoPid = argv[2] ;
		g_buzuoYmd0 = atof(argv[3]) ;
        g_buzuoYmd1 = atof(argv[4]) ;
        string extraCmd = argv[5] ;//no use
        spdlog::info("start for buzuo_plotonly from {} to {}.", g_buzuoYmd0,g_buzuoYmd1) ;
    }else if( argc!=2 ) {
        cout<<"unknow mode, params is "<<argc<<endl ;
        return 11 ;
    }


	GDALAllRegister();


	///////////////////////////////////////////////////////////////////////
	//
	//
	//  Initializing parameters.
	//
	///////////////////////////////////////////////////////////////////////

	g_pixelarea = 16;//4km*4km

    spdlog::info("parsing {}" , g_jsonfile) ;
    loadJsonParams(
            g_jsonfile, 
            g_gpp ,//global params
            g_fy4dayPP, 
            g_imsdayPP, 
            g_level3ppVec,
            g_jupingppVec,
            g_qhc
    ) ;
    g_gdaltranslate = g_gpp.gdaltranslate ;
    g_gdalwarp = g_gpp.gdalwarp ;
    g_plotcmd = g_gpp.plotcmd ;
    
    cout<<"parsing json done."<<endl ;
    
    
    //log
    //2023-1-10 comment 
    //auto dailylogger = spdlog::daily_logger_mt("l", g_gpp.log , 2, 0) ;//2:00 am

    //2023-1-10 use spdlog write into console only. log file use airflow.
    auto dailylogger = spdlog::stdout_color_mt("cl") ;
    spdlog::flush_every(std::chrono::seconds(5) ) ;
    spdlog::set_default_logger(dailylogger );

    //
    string mysqlConnError ;
    spdlog::info("mysql conn params:{} {} {} {}" , g_gpp.host, g_gpp.user, g_gpp.pwd, g_gpp.db);
    bool connOk = g_db.connect( g_gpp.host, g_gpp.user, g_gpp.pwd, g_gpp.db, mysqlConnError);
    if( connOk==false ){
        cout<<"mysql conn failed:"<<mysqlConnError<<endl ;
        spdlog::warn("mysql conn failed:{} ", mysqlConnError ) ;
    }


	cout<<"////////////////////////////////////////"<<endl ;
	cout<<"//"<<endl ;
	cout<<"//"<<endl ;
	cout<<"//  FY4A-SNC snow daily"<<endl ;
	cout<<"//"<<endl ;
	cout<<"////////////////////////////////////////"<<endl ;
    //记得要在87 上把fy4snc预处理程序吊起来，这边不再做预处理了
	if( g_runMode == RUNMODE_AUTO &&  g_fy4dayPP.enabled == 1 )
	{
		processFy4ImsDaily(g_fy4dayPP,"fy4");
	}
	else if( needBuzuo(g_buzuoPid , g_fy4dayPP) ==true )
    {
		processFy4ImsDaily(g_fy4dayPP,"fy4");        
    }else{
        cout<<"runmode:"<<g_runMode<<", fy4day.enable:"<<g_fy4dayPP.enabled<<endl ;
        cout<<"no need fy4 day"<<endl ;
    }



	cout<<"////////////////////////////////////////"<<endl ;
	cout<<"//"<<endl ;
	cout<<"//"<<endl ;
	cout<<"//  IMS snow daily"<<endl ;
	cout<<"//"<<endl ;
	cout<<"////////////////////////////////////////"<<endl ;
	if( g_runMode == RUNMODE_AUTO && g_imsdayPP.enabled==1 )
	{
		processFy4ImsDaily( g_imsdayPP , "ims");
	}else if( needBuzuo(g_buzuoPid , g_imsdayPP)==true )
	{
		processFy4ImsDaily( g_imsdayPP ,"ims" );
	}else{
        cout<<"no need ims day"<<endl ;
    }



	cout<<"////////////////////////////////////////"<<endl ;
	cout<<"//"<<endl ;
	cout<<"//"<<endl ;
	cout<<"//  process All Level3 "<<endl ;
	cout<<"//"<<endl ;
	cout<<"////////////////////////////////////////"<<endl ;
	processLevel3( g_level3ppVec );


	cout<<"////////////////////////////////////////"<<endl ;
	cout<<"//"<<endl ;
	cout<<"//"<<endl ;
	cout<<"//  process All JuPing products "<<endl ;
	cout<<"//"<<endl ;
	cout<<"////////////////////////////////////////"<<endl ;
    for(int ipp = 0 ; ipp < g_jupingppVec.size(); ++ ipp ){
        ProductParam& pp = g_jupingppVec[ipp] ;
        if( g_runMode == RUNMODE_AUTO && pp.enabled==1 )
        {
            processJupingProduct( pp , g_qhc);
        }else if( needBuzuo(g_buzuoPid , pp)==true )
        {
            processJupingProduct( pp , g_qhc );
        }else{
            cout<<"no need process "<< pp.pname <<endl ;
        }
    }


	std::cout << "All done." << std::endl;
	printCurrentTime();
	return 0;
}

